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In this paper I will describe some results that have been recently obtained in the study of 
random Euclidean matrices, i.e. matrices that are functions of random points in Euclidean space. 
In the case of translation invariant matrices one generically finds a phase transition between a 
phonon phase and a saddle phase. If we apply these considerations to the study of the Hessian of 
the Hamiltonian of the particles of a fluid, we find that this phonon-saddle transition corresponds 
to the dynamical phase transition in glasses, that has been studied in the framework of the mode 
coupling approximation. The Boson peak observed in glasses at low temperature is a remanent 
of this transition. We finally present some recent results obtained with a new approach where 
one deeply uses some hidden supersymmetric properties of the problem. 

1 Introduction 

In the last years many people have worked on the problem of analytically computing the properties 
of Euclidean random matrices [T]-|15|. The problem can be formulated as follows. 

We consider a set of N points (xi) that are randomly distributed with some given distribu- 
tion. Two extreme examples are: 

• The x's are random independent points with a flat probability distribution, with density p. 

• The cc's are one of the many minima of a given Hamiltonian. 

In the simplest case of the first example, given a function f(x), we consider the N x N matrix: 



The problem consists in computing the properties of the eigenvalues and of the eigenvectors of M. 
Of course, for finite N they will depend on the instance of the problem (i.e. the actually choice of 
the x's), however system to system fluctuations for intensive quantities (e.g. the spectral density 
g(z)) disappear when we consider the limit N — > oo at fixed particle density p. 

The problem is not new; it has been carefully studied in the case where the positions of the 
particles are restricted to be on a lattice [IB]. The case where the particles can stay in arbitrary 
positions, that is relevant in the study of fluids, has been less studied, although in the past many 
papers have been written on the argument [T]-|15|. These off-lattice models present some technical 
(and also physical) differences with the more studied on-lattice models. 

There are many possible physical motivations for studying these models, that may be applied 
to electronic levels in amorphous systems, very diluted impurities, spectrum of vibrations in glasses 
(my interest in the subject comes from this last field). 

I will concentrate in these lectures on the spectral density and on the Green functions, trying 
to obtaining both the qualitative features of these quantities and to present reasonable accurate 
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computations (when possible). This task is not trivial because there is no evident case that can be 
used as starting point for doing a perturbation expansion. Our construction can be considered as 
a form of a mean field theory (taking care of some corrections to it): more sophisticated problems, 
like localization or corrections to the naive critical exponents will be only marginally discussed in 
these notes. 

A certain number of variations to equation (JJ) are interesting. For example we could consider 
the case where we add a fluctuating term on the diagonal: 

M ilk = k,k f( X i ~ X i) ~ f( X i ~ X k) ■ ( 2 ) 

j 

This fluctuating diagonal term has been constructed in such a way that 

E M ^ = °- ( 3 ) 

k 

Therefore the diagonal and the off-diagonal matrix elements are correlated. 

In general it may be convenient to associate a quadratic form to the matrix M: the quadratic 
form il~L[fa\) is defined as: 

H[fa\ = ^J2fafaM i>k . (4) 

i,k 

In the first case we considered, eq. we have that: 

= H f( x i ~ x k)fa<Pk • (5) 
i,k 

In this second case, eq. (J2J), the associated quadratic form is given by 

H[fa = ]T fafaM i>k = ~J2 f( x i - - <Pk? ■ (6) 

Here the matrix M is non-negative if the function / is non-negative. The matrix M has always a zero 
eigenvalue as consequence of the invariance of the quadratic form under the symmetry fa — > fa + A; 
the presence of this symmetry has deep consequences on the properties of the spectrum of M and, 
as we shall see later, phonons are present. Many of the tricky points in the analytic study are 
connected to the preserving and using this symmetry in the computations. 

In the same spirit we can consider a two-body potential V(x) and we can introduce the Hamil- 
tonian 

H[x] = l -Y. V ^- x k) ■ (7) 

i,k 

We can construct the 3iV x 3N Hessian matrix 

d 2 H 

M itk = = S i>k J2 V"{ Xi -x 3 )- V"( Xl -x k ) , (8) 

dxix k ^ 

where for simplicity we have not indicated space indices. Also here we are interested in the com- 
putation of the spectrum of M . The translational invariance of the Hamiltonian implies that the 
associated quadratic form is invariant under the symmetry fa — > fa + A and a phonon-like behavior 
may be expected. 

This tensorial case, especially when the distribution of the x's is related to the potential V, is 
the most interesting from the physical point of view (especially for its mobility edge |16[ I15j ) . Here 
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we stick to the much simpler question of the computation of the spectrum of M. We shall develop 
a field theory for this problem, check it at high and low densities, and use a Hartree type method. 

Our aim is to get both a qualitative understanding of the main properties of the spectrum and 
of the eigenvalues and a quantitate, as accurate as possible, analytic evaluation of these properties. 
Quantitative accurate results are also needed because the computation of the spectral density is a 
crucial step in the microscopic computation of the thermodynamic and of the dynamic properties 
of glass forming systems. In some sense this approach can be considered as an alternative route to 
obtain mode-coupling like results |17) , the main difference being that in the mode coupling approach 
one uses a coarse grained approach and the hydro dynamical equations, while here we take a fully 
microscopic point of view. 

Many variations can be done on the distribution of points, each one has its distinctive flavor: 

• The distribution of the points x is flat and the points are uncorrelated: they are uniformly 
distributed in a cube of size L and their number is N = pL 3 . Here, as in the following cases, 
we are interested in the thermodynamic limit where both L and N go to infinity at fixed p. 

• The points x are distributed with a distribution that is proportional to exp(—/3H[x]) where 
H[x] is a given function. 

• The points x are one of the many solutions of the equation dH/dxi = 0. This last problem 
may be generalized by weighting each stationary point of H with an appropriate weight. 

The last two cases are particularly interesting when 



and consequently 



If this happens the distribution of points and the matrix are related and the theoretical challenge 
is to use this relation in an effective way. 

In the second section of these lectures, after the introduction we will present the basic definition 
of the spectrum, Green functions, structure functions and related quantities. In the third section 
we will discuss the physical motivation for this study and we will present a first introduction to 
the replica method. In the fourth section we will give a few example how a field theory can be 
constructed in such way that it describes the properties of randomly distributed points. In the 
fifth section we will discuss in details the simples possible model for random Euclidean matrices, 
presenting both the analytic approach and some numerical simulations. In the sixth section we shall 
present a similar analysis in a more complex and more physically relevant case, where phonons are 
present due to translational invariance. Finally, in the last section we present some recent results 
that have been obtained in the case of correlated points. 

2 Basic definitions 

The basic object that we will calculate are the matrix element of the resolvent 
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(we use the notation M[x\ in order to stress that the matrix M depends on all the points [x\. 
Sometimes we will suppress the specification "[x]", where it is obvious. We can define the sample 
dependent Green function 



1,3 



G[x}( yi ,y 2 ,z) = ^25(y 1 - Xi)8{y 2 - Xj)G[x}(z) id 

1 

z - M[x] 



^<%i - Xi)S(y 2 - xj) 



1,3 



(12) 



The quantity that we are interested to compute analytically are the sample averages of the Green 
function, i.e. 

G(yi,y 2 ,z) = G[x]( yi ,y 2 ,z) , (13) 

where the overline denotes the average over the random position of the points. 

If the problem is translational invariant, after the average on the positions of the points [x] we 
have 

G{yi,y 2 ,z) =G(y 1 -y 2 ,z) , (14) 

where the function G(x, z) is smooth apart from a delta function at x = 0. It is convenient to 
consider the Fourier transform of G(x, z), i.e. 



G{p,z) = — Y^jv&i-*,) 



1 



z-M 



dxdye^^G^ix.y.z) = J dxe ipx G(x, z) , 



(15) 



The computation of the function G(p, z) will one of the main goals of these lectures. It is con- 
venient to introduce the so called dynamical structure factor (that in many case can be observed 
experimentally) that here is defined as l : 



S e {p,E) = lim G(p,E + irj) . 

7T ?7^0+ 



(16) 



One must remember that the structure function is usually defined using as variable u> = \[E: 

S(p,u) =2cu S e (p,uj 2 ), (17) 
The resolvent will also give us access to the density of states of the system: 

1 



N r 



9e(E) 



lim 

' i=i 

1 



1 



E + irj- M 



= lim lim G(p, E + iri) . 

7T r^0+ 



(18) 



3 Physical motivations 

There are many physical motivations for studying these problems, apart from the purely mathe- 
matical interest. Of course different applications will lead to different forms of the problems. In 
this notes I will concentrate my attention on the simplest models, skipping the complications (that 
are not well understood) of more complex and interesting models. 

1 We shall often make use of the distribution identity l/(x + i0 + ) = P(l/x) — br6(x) where P(l/x) denotes the 
principal part, i.e. l/(x + i0 + ) + l/(x + iCT) = 2P(l/x). 
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3.1 Impurities 

We can consider a model where there are some impurities in a solid that are localized at the points 
x's. There are many physical implementation of the same model; here we only present two cases: 

• There may be electrons on the impurities and the amplitude for hopping from an impurity (i) 
to an other impurity (k) is proportional to f(xi — xj-)- The electron density is low, so that 
the electron-electron interaction can be neglected. 

• There is a local variable (pi associated to the impurity (e.g. the magnetization) and the effective 
Hamiltonian at small magnetizations is given by eq.(jSJ). 

In both cases it is clear that the physical behavior is related to the properties of the matrix M 
that is of the form discussed in the introduction. If the concentration of the impurities is small 
(and the impurities are embedded in an amorphous system) the positions of the impurities may be 
considered randomly distributed in space. 



3.2 Random walks 

Let us assume that there are random points in the space and that there is a particle that hops from 
one point to the other with probability per unit time given by f(xi — Xf.). 

The evolution equation for the probabilities of finding a particle at the point i is given by 

dP 

k k k 

where the matrix M is given by eq. (J2J). 

Let us call P(x, t) the probability that a particle starting from a random point at time arrives 
at at time t. It is evident that after averaging of the different systems (or in a large system after 
averaging over the random starting points) we have that 

J dtexp(—tz) J dx exp(ipx)P(x, t) = — G(p, — z) , (20) 

so that the function P(x, t) can be reconstructed from the knowledge of the function G(p, z) defined 
in eq. (fT3)l . 



3.3 Instantaneous modes 

We suppose that the points x's are distributed according to the probability distribution 

exp(-PH[x]) , (21) 

where H[x) is for example a two body interaction given by eq.(J7J) and (3 = l/(kT). 

The eigenvectors of the Hessian of the Hamiltonian (see eq.©) are called instantaneous normal 
modes (INN). The behavior of the INN at low temperature, especially for systems that have a glass 
transition |20l I19j . is very interesting |21| . In particular there have been many conjectures that 
relate the properties of the spectrum of INN to the dynamic of the system when approaching the 
glass transition. 

However it has been realized in these years that much more interesting quantities are the instan- 
taneous saddle modes (ISN) [221 12H1 12E1 EU • Given an equilibrium configuration x e at temperature 
T one defined the inherent saddle x s corresponding to that configuration x e in the following way. 
We consider all the saddles of the Hamiltonian, i.e. all the solutions of the equations 
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The inherent saddle x s is the nearest saddle to the configuration x e . The Boltzmann distribution at 
temperature T induces a probability distribution on the space of all saddles. The physical interesting 
quantities are the properties of the Hessian around these saddles. It turns out that they have a very 
interesting behavior near the glass transition in the case of glass forming materials. 

It is clear that the computation of the spectrum in these case is much more complex: we have 
both to be able to compute the correlations among the points in this non-trivial case and to take 
care of the (often only partially known) correlations. 

A possible variation on the same theme consists in consider the ensemble of all saddles of given 
energy E = H[x\: there are indications that also in this case there may be a phase transition in 
the properties of the eigenvalues and eigenvectors when E changes and this phase transition may 
be correlated to some experimentally measured properties of glassy systems, i.e. the Boson Peak, 
as we will see later. 

4 Field theory 

4.1 Replicated Field Theory 

The resolvent G(p, z) defined in eq. (|15|) can be written as a propagator of a Gaussian theory: 

z[j] = |n^ ex p(-^E^^- M )i m ^+E^ j o > 

k \ Im I / 

G ^ z ) = ^ErTTT ex P(^-^))ln^[J] | j=Q , (23) 

ij 1 3 

where the over line denotes the average over the distribution P[x] of the positions of the particles. 
We have been sloppy in writing the Gaussian integrals (that in general are not convergent): one 
should suppose that the (f) integrals go from — aoo to +aoo where a 2 = —i and z has a positive 
(albeit infinitesimal) imaginary part in order to make the integrals convergent |16l I18j . This choice 
of the contour is crucial for obtaining the non compact symmetry group for localization |18| . but it 
is not important for the density of the states. 

In other words, neglecting the imaginary factors, we can introduce a probability distribution 
over the field fa, 

P[4>] oc J] dfa exp ( ~ <h ( z ~ M \m <Pm) ■ (24) 

k \ Im / 

If we denote by (•) the average with respect to this probability (please notice that here everything 
depends on [x] also if this dependence is not explicitly written) we obtain the simple result: 

= G hk . (25) 

A problem with this approach is that the normalization factor of the probability depends on [x] 
and it is not simple to get the x-averages of the expectation values. Replicas are very useful in this 
case. They are an extremely efficient way to put the dirty under a simple carpet (on could also 
use the supersymmetric approach of Efetov where one puts the same dirty under a more complex 
carpet [To]^. 

For example the logarithm in eq. (|23|) is best dealt with using the replica trick [Tll2ll25]: 

In Z[J} = lim-(Z n [J]-l) . (26) 
?wo n 

The resolvent can then be computed from the n-th power of Z, that can be written using n replicas, 
cf)^ (a = 1, 2, . . . ,n) of the Gaussian variables of ea.(|23|l: 
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G(p,z) = lim GW(p,z) , 

NG {n \p,z)= (27) 



Indeed 




and the physically interesting case is obtained only when n = 0. 

In this way one obtains a 0(n) symmetric field theory. It well known that an 0(n) symmetric 
theory is equivalent to an other theory invariant under the 0(n + 2|1) symmetry (PHI)) where this 
group is defined as the one that leaves invariant the quantity 

£ (<n 2 + ^, (29) 

a=l,n+2 

where ip is a Fermionic field. Those who hate analytic continuations, can us the group 0(2|1) at 
the place of the more fancy O(0) that is defined only as the analytic continuation to n = of the 
more familiar 0(n) groups. 

Up to this point everything is quite general. We still have to do the average over the random 
positions. This can be done in an explicit or compact way in some of the previous mentioned cases 
and one obtains a field theory that can be studied using the appropriate approximations. Before 
doing this it is may be convenient to recall a simpler result, i.e. how to write a field theory for the 
partition function of a fluid. 



4.2 The partition function of a fluid 

Let us consider a system with N (variable) particles characterized by a classical Hamiltonian Hn[x] 
where the variable x denote the ./V positions. In the simplest case the particles can move only in a 
finite dimensional region (a box B of volume Vb) and they have only two body interactions: 

Hn[x] = ^V(xi-x k ) . (30) 

i,k 

At given (3 the canonical partition function can be written as 

Qn = / dxi . . . dxjs[ exp(— j3Hjs[[x\) , (31) 
Jb 

while the gran-canonical partition function is given by 

2(0 = £^QJV, (32) 

N 

where Q is the fugacity. 
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We aim now to find out a field theoretical representation of the Q(C)- This can be formally done 
:28_ (neglecting the problems related to the convergence of Gaussian integrals) by writing 

exp (^-^^2V{xi- x k ) 

N~ l J daexp(^ J dxdyV~ l {x - y)a{x)a{y) + ^<r(x i ) S j , (33) 

where M is an appropriate normalization factor such that 

M~ l J da exp J dxdyV~ 1 {x - y)a(x)a(y) + j J{x)a(x)^j = 

exp (- J dxdyPV(x - y)J{x)J{y^j . (34) 
The reader should notice that V~ 1 (x — y) is formally defined by the relation 

dzV- 1 {x-z)V(z-y) = 5(x-y) . (35) 

The relation eq. ()33|) trivially follows from the previous equations if we put 

J(x) = Y^6(x-x i ) . (36) 

i 

It is convenient to us the notation 

dn[a] = ATWexp (-L J dxdyV~\x - y)a{x)a(y)^j . (37) 
With this notation we find 

Qn = J dp[a] (j dxexp(a(x)^j . (38) 

We finally get 

Q(C) = J dp[a]exp (( J dxexp(a(x))) . (39) 

People addict with field theory can use this compact representation in order to derive the virial 
expansion, i.e. the expansion of the partition in powers of the fugacity at fixed j3. 

The same result can be obtained in a more straightforward way |29| by using a functional integral 
representation for the delta function: 

exp V ( Xi ~ Xk ) 

j d[p]5 F p(x) - S(x-Xi) expf-^ j dxdy p(x) p(y)V \x - y) \ , (40) 
J [ i=l,JV J v A J y 

where 8f stands for a functional Dirac delta: 

5 F [f] = J da exp (i J dxf{x)\(x)\ . (41) 
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We thus find 



Q N = J dxidx N J d[p]d[\] (42) 
exp [ — ^ / dxdy p{x) p{y)V {x — y) + i\(x) p{x) — i ^ \{x 

\ J i=l,N 



d[p]d[\] exp ( - ^ / dxdyp{x)p(y)V{x - y) + iX(x)p(x) 



2 



. N 

dx exp(iA(x) 



At the end of the day we get 

Q(Q = J d[X]exp J dxdyX(x)X(y)V~ 1 (x - y) + ( J dxexp(i\(x)\ . (43) 

where, by doing the Gaussian integral over the p, we have recover the previous formula ea. (|33jl . if 

we set iX(x) = a{x). 

In order to compute the density and its correlations it is convenient to introduce an external 
field 

J2U( Xi ). (44) 

i 

This field is useful because its derivatives of the partition function give the density correlations. As 
before the partition function is given by 

Q(C\U) = J dp[a]exp (( J dx exp(a(x) + f3U(x))J . (45) 

In this way one finds that the density and the two particle correlations are given by: 

p=C£5(x-x l )) = ((exp(a(0))) 

i 

C(x) = 5( Xl - 0)5(x k - x)} = C 2 (exp(a(0) + a(x))) . (46) 

i,k 

In the low density limit (i.e. C near to zero) one finds 

p = Cexp(-/?y(0)), 
C{x) = p 2 exp(-(3V(x)) . (47) 

that is the starting point of the virial expansion. 

5 The simplest case 

In this section we shall be concerned with simple case for random Euclidean matrices, where the 
positions of the particles are completely random, chosen with a probability P{x) = p/Vb- 
We will consider here the simplest case, where 

M t , k = f( Xi - x k ) (48) 

and fix) is bounded, fast decreasing function at infinity (e.g. exp(— x 2 /2)). 

We shall study the field-theory perturbatively in the inverse density of particles, 1/p. The 
zexo th order of this expansion (to be calculated in subsection l5.1|) corresponds to the limit of infinite 
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density, where the system is equivalent to an elastic medium. In this limit the resolvent (neglecting 
p independent terms) is extremely simple: 

G(p,z) = (49) 

In the above expression f(p) is the Fourier transform of the function /, that due to its spherical 
symmetry, is a function of only (p 2 ) . We see that the dynamical structure function has two delta 



functions at frequencies u> = ±y pf(p). 

It is then clear that Eq. (|49|) represents the undamped propagation of plane waves in our harmonic 
medium, with a dispersion relation controlled by the function / . The high order corrections to 
Eq. (j49|) (that vanishes as 1/p) will be calculated later. They take the form of a complex self-energy, 
£(p,3), yielding 

G( P ,z) = rY i— r, (50) 

E{ ' K(E-e{p)-Re^{p,E)) 2 + (ImZ(p, E)) 2 ' 1 ' 

The dynamical structure factor, is no longer a delta function, but close to its maxima it has a 
Lorentzian shape. From eq. (|51|) we see that the real part of the self-energy renormalizes the 
dispersion relation. The width of the spectral line is instead controlled by the imaginary part of the 
self-energy. 

subsectionField theory Let us first consider the spectrum: it can be computed from the imaginary 
part of the trace of the resolvent: 

R ^ = ^jhi ■ < 52 > 

where the overline denotes the average over the positions X{. 

It is possible to compute the resolvent from a field theory written using a replica approach. We 
shall compute H^v = det(z — M)~ n / 2 , and deduce from it the resolvent by using the replica limit 
n -> 0. 

It is easy to show that one can write S^v as a partition function over replicated fields (fif, where 
i € {1...N}, a e {l...n}: 

N , „ N n 

ax; 



exp 



n^/nn^ 

i=l J i=l a=l 

(-fEc^ + sE/^-^w) ■ (53) 

\ i,a i,j,a / 

In order to simplify the previous equations let us introduce the Bosonic 2 fields ^ a (x) = YliLi 4% S(x— 
xi) together with their respective Lagrange multiplier fields ip a (x). More precisely we introduce a 
functional delta function: 

Sf (^a(x) - 4>tK x ~ x i)^j = 
J d[ip a ] exp (i J dx^ a {x) (ip a {x) - E ~ ) ■ ( 54 ) 



2 In spite of a long tradition, here the fields ip are Bosons, not Fermions. Fermionic fields will appear only in the 
last section of these note. 
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One can integrate out the 4> variables, leading to the following field theory for Sjv, where we have 
neglected all the z independent factors that disappear when n = 0: 



J D[ip a ,ip a }A N exp(S ) (55) 



where 



S o = *H/ dx Mx)ifi a (x) + J dxdy i; a {x)f(x -y)i; a (y) , 



A = exp 



(56) 



It is convenient to go to a grand canonical formulation for the disorder: we consider an ensemble 
of samples with varying number of (N), and compute the grand canonical partition function Z(C) = 
T,n=o z>n( N /N\ that is equal to: 

Z = f D[^ a ] exp (So + (A) ; . (57) 

In the n — > limit, one finds that p = Q . The definition of the field tp implies that the correlation 
is given (at x / 0) by: 

G(x, z) = J2 S(x - Xi)8{ Xj ) (—^f) = (Mx)MO)) (58) 



A simple computation (taking care also of the contribution in x = 0) gives 

G(p, z) = - lim -J— f d d xd d y exp(ip(x - y)) $W (x)^ (y)) . (59) 

so that the average Green function can be recovered from the knowledge of the tp propagator. 
Notice that we can also integrate out the ip field thus replacing So by S' , where 

S ° = \llJ dxd y ^ a ( x )f' l {x - yW(y), (60) 

where / _1 is the integral operator that is the inverse of /. 

The expression (|57|) is our basic field theory representation. We shall denote by brackets the 
expectation value of any observable with the action So + Si. As usual with the replica method we 
have traded the disorder for an interacting replicated system. 

It may be convenient to present a different derivation of the previous result: we write 

exp ( - Z - J>«) 2 + I £ f( Xi - ZjWfl 

/ d[u] exp (-- Y^f dxdy u(x) a u(y) a f-\x - y) + g (^(x^x^ - ^) 2 ) j • (61) 
If we collect all the terms at a given site i, we get 

J Hdtfexp " § (€) 2 )) • (62) 
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The 4> integrals are Gaussian and they can be done: for n = one remains with 



exp ]^2u(xi)l/zJ (63) 

and one recovers easily the previous result ea. Q56JI . 

The basic properties of the field theory are related to the properties of the original problem in 
a straightforward way. We have seen that the average number of particles is related to £ through 
N = (V{A), so that one gets £ = p in the n — > limit because (A) = 1. From the generalized 
partition function Z, one can get the resolvent R(z) through: 

R(z) = - lim 4^ . (64) 

5.1 High density expansion 

Let us first show how this field theory can be used to derive a high density expansion. 

At this end it is convenient to rescale z as z = pz and to study the limit p — > oo at fixed z. 
Neglecting addictive constants, the action, S, can be expanded as: 

i £ / dxdy r(x)f- 1 (x - yW{y) - 1 J dxj^ M^f + 0(l/(pz 2 ) , (65) 

Gathering the various contributions and taking care of the addictive constants, one gets for the 
resolvent: 

11/"../ 1 1 



pR(z) = - + - dkl = . (66) 

z pJ \(z-f(k)) z) 

One can study with this method the eigenvalue density for eigenvalues |A| ~ O(p), by setting 
z = A + ie and computing the imaginary part of the resolvent in the small e limit. For p — > oo the 
leading term gives a trivial result for the eigenvalue density <?(A) i.e. 

g(X) = S(X) . (67) 

Including the leading large p correction that we have just computed, we find that g(X) develops, 
away from the peak at A ~ 0, a component of the form: 

g(X) ~ - / dk8(X - f(k)) . (68) 
pJ 

The continuos part of the spectrum can also be derived from the following simple argument. We 
suppose that the eigenvalue U{ is a smooth function of Xj. We can thus write: 

fix - x k )u){x k ) w p I dyfix - y)uj(y) = pXu;(x) (69) 

k J 

and the eigenvalues are the same of the integral operator with kernel /. 

This argument holds if the discrete sum in eq. (|69|) samples correctly the continuous integral. 
This will be the case only when the density p is large enough that the function u>(x) doesn't oscillate 
too much from one point Xj to a neighboring one. This condition in momentum space imposes that 
the spatial frequency |fc| be small enough: \k\ <C p l ^ d ( p l ^ d is the inverse of the average interparticle 
distance) . 

The same condition is present in the field theory derivation. We assume that f(k) decreases at 
large k, and we call ku the typical range of k below which f(k) can be neglected. Let us consider 
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the corrections of order p _1 in eq. Q66[) . It is clear that, provided z is away from 0, the ratio of the 
correction term to the leading one is of order kf^p -1 , and the condition that the correction be small 
is just identical to the previous one. The large density corrections near to the peak z = cannot 
be studied with this method. 

In conclusion the large p expansion gives reasonable results at non-zero z but it does not produce 
a well behaved spectrum around z = 0. In the next sections we shall see how to do it. 

The reader may wonder if there is a low density expansion: the answer is positive however 
it is more complex than the high density expansion and it cannot be discussed in these notes for 
lack of space. 

5.2 A direct approach 

Let us compute directly what happens in the region where p is very large. We first make a very 
general observation: if the points x are random and uncorrelated we have that 

J2 A (xi) = P J dxA(x) , 

i 

^B(xi,x k )= B(xi,x k ) + ^2B(xi,Xi) = 

i,k i,k\i^k i 

p 2 J dxdyB(x,y) + p J dxB(x,x) (70) 

When the density is very high the second term can be neglected and we can approximate multiple 
sum with multiple integrals, neglecting the contribution coming from coinciding points |2|IH1IZ|- 

We can apply these ideas to the high density expansion of the Green function in momentum 
space. Using a simple Taylor expansion in 1/z we can write: 

G{p,z) = -Y J (-z)- R MHp) (71) 

Z R 

where for example 

M 3 (p) = N" 1 Y f( x k ~ x h)f(x kl - x k2 )f(x k2 - x ks ) exp(ip(x fco - x ks ) (72) 

k ,kik2,k 3 

The contribution of all different points gives the large p limit, and taking care of the contribution 
of pairs of points that are not equal gives the subleading corrections. In principle everything can 
be computed with this kind of approach that is more explicit than the field theory, although the 
combinatorics may become rather difficult for complex diagrams and the number of contributions 
become quite large when one consider higher order corrections [SI [7]. 

Generally speaking the leading contribution involve only convolutions and is trivial in momentum 
space, while at the order l/rho m , there are m loops for diagrams that have m intersections: the 
high p expansion is also a loop expansion. 

5.3 A variational approximation 

In order to elaborate a general approximation for the spectrum, that should produce sensible results 
for the whole spectrum, we can use a standard Gaussian variational approximation in the field theory 
representation, that, depending on the context, appears under various names in the literature, like 
Hartree-Fock approximation the Random Phase Approximation or CPA. 

Here we show the implementation in this particular case. By changing ip — > itp in the represen- 
tation for the partition function, we obtain: 

Z = J D[fr]exp(S)) , (73) 



13 



where 



s = -\J dxd y Y,4> a (x)rHx,yW(y) 

+pz- n / 2 J dxexp (^E^O*) 2 ) • ( 74 ) 
We look for the best quadratic action 

S v = -(1/2) J2 J dxdyG-J(x,yW{x)^ h (y) (75) 

ab 

that approximates the full interacting problem. 

This procedure is well known in the literature and it gives the following result. If we have a field 
theory with only one field and with action 



dpD(p)(f){p) 2 + / dxV(<f>(x)), (76) 



the propagator is given by 



where 



0(rt " dwTe ' (77) 



S = . (78) 



In other words £ is a momentum independent self-energy that is equal to the expectation value of 
V"(4>) in theory where the field (ft is Gaussian and has a propagator equal to G. 

In the present case, there are some extra complications due to the presence of indices; the 
appropriate form of the propagator is G a b(p) = 5 a bG(p). After same easy computations, one finds 
that G(p) satisfies the self consistency equation: 

G(p)= f _ l( \ y , S = -^- 7 — (79) 

/ \P) ~ £ z- J dkG(k) 

and the resolvent is given by: 

Riz) = ^ . (80) 

V ' z-JdkG(k) V ' 

Formulas (|79l80j) provide a closed set of equations that allow us to compute the Gaussian vari- 
ational approximation to the spectrum for any values of z and the density; if the integral would 
be dominated by a single momentum, we would recover the usual Dyson semi-circle distribution 
for fully random matrices . In sect. 15.41 we shall compare this approximation to some numerical 
estimates of the spectrum for various functions / and densities. The variational approximation 
correspond to the computation of the sum of the so called tadpole diagrams, that can be often done 
a compact way. 

The main advantage of the variational approximation is that the singularity of the spectrum at 
z = is removed and the spectrum starts to have a more reasonable shape. 

Another partial resummation of the p expansion can also be done in the following way: if one 
neglects the triangular-like correlations between the distances of the points (an approximation that 
becomes correct in large dimensions, provided the function / is rescaled properly with dimension), 
the problem maps onto that of a diluted random matrix with independent elements. This problem 
can be studied explicitly using the methods of |3U1 131j . It leads to integral equations that can be 
solved numerically. The equations one gets correspond to the first order in the virial expansion, 
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Figure 1: Density of eigenvalues of a Euclidean Random Matrix in three dimensions, density p = 1. 
The function / is f(x) = (27r)~ 3//2 exp(— x 2 /2), and the matrix is defined from eq.(|48|). The full 
line is the result of a numerical simulation with N = 800 points, averaged over 100 samples. The 
dashed line is the result from the high density expansion. The dash-dotted line is the result from 
the Gaussian variational approximation (RPA) to the field theory. 

where one introduces as variational parameter the local probability distribution of the field <j> |32j . 
The discussion of this point is interesting and it is connected to the low-density expansion: however 
it cannot be done here. 

An other very interesting point that we will ignore is the computation of the tails in the prob- 
ability distribution of the spectrum representation leading to compact expressions for the behavior 
of the tail [S3]. 

5.4 Some numerical experiments 

In this section we present the result of some numerical experiments and the comparison with the 
theoretical results (we also slightly change the notation ,i.e. we set z = A + ie. 

We first remark that for a function f{x) that has a positive Fourier transform f(k), the spectrum 
is concentrated an the positive axis. Indeed, if we call uJi a normalized eigenvector of the matrix M 
defined in (|48|1. with eigenvalue A, one has: 

Uif(xi - Xj)ujj = A , (81) 

ij 

and the positivity of the Fourier transform of / implies that A > 0. 

We have studied numerically the problem in dimension d = 3 with the Gaussian function f(x) = 
(27r) -3 / 2 exp(— x 2 /2). In this Gaussian case the high density approximation gives a spectrum 

5(A) ~ -\ \ Q log Q V2 6(p - A) + C5(X) (82) 

Notice that this spectrum is supposed to hold away from the small A peak, and in fact it is not 
normalizable at small A. 
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Figure 2: Density of the logarithm (in base 10) of the eigenvalues of a Euclidean Random Matrix 
in three dimensions, density p = 1 (same data of the previous figure). The dashed line is the 
result from the high density expansion. The dotted line is the result from the Gaussian variational 
approximation (RPA) to the field theory. 



It is possible to do a variational approximation computation taking care of all the non-linear 
terms in the action for the fields. This corresponds (as usually) a a resummation of a selected class 
of diagrams. After some computations, one finds that one needs to solve, given z = A — ie, the 
following equations for C(z) = a(z) + ib(z): 



a 1 r°° 2 e k > z - u 

P 0? + b 2 + 2^2 J Q ( e fc2/ 2 _ a ) 2 + 6 2 

b b f°° 2 1 

VT^^X fcc V/ 2 - a ) 2 + ^ (83) 

One needs to find a solution in the limit where e — > 0. 

In fig. (|1I2I) . we plot the obtained spectrum, averaged over 100 realizations, for N = 800 points 
at density p = 1 (We checked that with a different number of points the spectrum is similar). Also 
shown are the high density approximation 1)82(1 . and the result from the variational approximation. 
We see from fig.Q that the part of the spectrum A G [0.2, 1.5] is rather well reproduced from both 
approximations, although the variational method does a better job at matching the upper edge. 
On the other hand the probability distribution of the logarithm of the eigenvalues (figEJ) makes it 
clear that the high density approximation is not valid at small eigenvalues, while the variational 
approximation gives a sensible result. One drawback of the variational approximation, though, is 
that it always produces sharp bands with a square root singularity, in contrast to the tails that are 
seen numerically. 

In figEl we plot the obtained spectrum, averaged over 200 realizations, for N = 800 points at 
density p = 0.1. We show also a low density approximation (that we do not describe here J2J), and 
the result from the variational approximation. We see from fig.® that this value of p = 0.1 is 
more in the low density regime, and in particular there exists a peak around A = /(0) due to the 
isolated clusters containing small number of points. The variational approximation gives the main 
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Figure 3: Density of eigenvalues of a Euclidean Random Matrix in three dimensions, density 
p = 0.1. The function / is f(x) = (27r)~ 3 / 2 exp(— x 2 /2), and the matrix is defined from eq. (|48|) 
with u = 0. The full line is the result of a numerical simulation with N = 800 points, averaged 
over 200 samples. The dashed line is the result from a low density expansion introduced in |2j. 
The dash-dotted line is the result from the Gaussian variational approximation (RPA) to the field 
theory. 

orders of magnitude of the distribution, but it is not able to reproduce the details of the spectrum, 
in particular the peak due to small clusters. On the other hand the leading term of a low density 
approximation (introduced in [2j) gives a poor approximation the the overall form of the spectrum. 
One should use an approach where the advantages of both methods are combined together. 

6 Phonons 

6.1 Physical Motivations 

Inelastic X-ray scattering (IXS) experiments and inelastic neutron scattering on structural glasses 
and supercooled liquids provided useful information on the dynamics of their amorphous structure, 
at frequencies larger than 0.1 THz (see for example |34j and references therein). Those experiments 
show a regime, when the wavelength of the plane wave is comparable with the inter-particle distance, 
where the vibrational spectrum can be understood in terms of propagation of quasi-elastic sound 
waves, the so-called high frequency sound. This high-frequency sound has also been observed in 
molecular dynamical simulations of strong and fragile liquids, and it displays several rather universal 
features. In particular, a peak is observed in the dynamical structure factor at a frequency that 
depends linearly on the exchanged momentum p, in the region O.lpo — 1-Opoj Po being the position 
of the first maximum in the static structure factor. When extrapolated to zero momentum, this 
linear dispersion relation yields the macroscopic speed of sound. The width of the spectral line, T 
is well fitted by 

T(p)=Ap e , x«2, (84) 

with A displaying a very mild (if any) temperature dependence. Moreover, the same scaling of T 
has been found in harmonic Lenhard- Jones glasses [HS]) an d one can safely conclude that the p 2 
broadening of the high-frequency sound is due to small oscillations in the harmonic approximation. 
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In these context other interesting problems related with the high-frequency vibrational excitations 
of these topologically disordered systems ^5] regard the origin of the Boson peak or the importance 
of localization properties to understand the dynamics of supercooled liquids [37j . 

The variety of materials where the p 2 broadening appears suggests a straightforward physical 
motivation. However, the simplest conceivable approximation, a wave propagating on an elastic 
medium in the presence of random scatterers, yields Rayleigh dispersion: T oc p 4 . This result is 
very robust: as soon as one assumes the presence of an underlying medium where the sound waves 
would propagate undisturbed, as in the disordered-solid model (HHl OBI EHJ , the p 4 scaling appears 
even if one studies the interaction with the scatterers non-perturbatively jlU]. When the distinction 
between the propagating medium and the scatterers is meaningless (as it happens for topologically 
disordered systems), the p 2 scaling is recovered. 

We want to investigate the problem from the point of view of statistical mechanics of random 
matrices, by assuming that vibrations are the only motions allowed in the system. The formalism 
we shall introduce, however, is not limited to the investigation of the high frequency sound and it 
could be straightforwardly applied in different physical contexts. 

Let us look more carefully at the relation between vibrational dynamics in glasses and random 
matrices. The dynamical structure factor for a system of ./V identical particles is defined as: 

S(p, u) = ^Y.j dt eiUt ( e ip - (rj(t) - ri(0)) ) , (85) 

where (. . .) denotes the average over the particles positions r,- with the canonical ensemble. 

Here it is convenient to consider the normal modes of the glass or the supercooled liquid, usually 
called instantaneous normal modes (INM) because they are supposed to describe the short time 
dynamics of the particles in the liquid phase [25- O ne studies the displacements u around the 
random positions x, by writing the position of the i-th particle as ri(t) = re, + Ui(t), and linearizing 
the equations of motion. Then one is naturally lead to consider the spectrum of the Hessian matrix 
of the potential, evaluated on the instantaneous position of the particles. Calling uj 2 the eigenvalues 
of the Hessian matrix and e n (i) the corresponding eigenvectors, the one excitation approximation 
to the S(p,uj) at non zero frequency is given in the classical limit by: 

k T N 

S^ l \p,uj) = -^Y,Qn(p)Hu-uj n ), (86) 
muj z — ' 

n=l 

Qui?) = p ■ e n (i) exp(ip ■ Xi)\ 2 . (87) 

i 

However one cannot always assume that all the normal modes have positive eigenvalues, negative 
eigenvalues representing the situation where some particles are moving away from the position x. 
Indeed, it has been suggested [37| that diffusion properties in supercooled liquids can be studied 
considering the localization properties of the normal modes of negative eigenvalues. 

A related and better defined problem is the study of normal modes at zero temperature, where 
the displacements u are taken around the rest positions. By assuming that this structure corresponds 
to one minimum of the potential energy, one can introduce a harmonic approximation where only 
the vibrations around these minima are considered, and all the dynamical information is encoded 
in the spectral properties of the Hessian matrix on the rest positions. The Hessian is in a first 
approximation a random matrix if these rest positions correspond to a glass phase. It has been shown 
using molecular dynamics simulation that below the experimental glass transition temperature 
the thermodynamical properties of typical strong glasses are in a good agreement with such an 
assumption. 

Therefore, the problem of the high-frequency dynamics of the system can be reduced, in its 
simplest version, to the consideration of random Euclidean matrices, where the entries are deter- 
ministic functions (the derivatives of the potential) of the random positions of the particles. As 
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far as the system has momentum conservation in our case, due to translational invariance, all the 
realizations of the random matrix have a three common normal mode with zero eigenvalue: the 
uniform translation of the system. 

An Euclidean matrix is determined by the particular deterministic function that we are consid- 
ering, and by the probabilistic distribution function of the particles, that in the INM case is given 
by the Boltzmann factor. However, for the sake of simplicity we shall concentrate here on the sim- 
plest kind of euclidean matrices without spatial correlations and we will neglect the vector indices 
of the displacement. We consider N particles placed at positions Xi,i = 1,2, ...,N inside a box, 
where periodic boundary conditions are applied. Then, the scalar euclidean matrices are given by 
eq.©, where f(x) is a scalar function depending on the distance between pairs of particles, and the 
positions {x} of the particles are drawn with a fiat probability distribution. Notice that the matrix 
(j2J) preserves translation invariance, since the uniform vector eo(i) = const is an eigenvector of M 
with zero eigenvalue. Since there are not internal indices (the particle displacements are restricted 
to be all collinear), we cannot separate longitudinal and transversal excitations. 

The dynamical structure factor for a scalar Euclidean matrix is given by 



S E (p,E) = Y,Q n (p)6(E - E n ) , 



(88) 



Qn{p) 

S(p,w) 



1 

N 



N 



e n (i)exp(ip 



i=l 



2u)Se{p, 



(89) 
(90) 



where the overline stands for the average over the particles position and we have given the definition 
either in the eigenvalue space (Se(p,E)) and in the frequency space (S(p,uj)). 



6.2 A more complex field theory representation 

The basic field theory representation is similar to the one of the previous section. The main com- 
plication is due to the presence of a diagonal term in the matrix M. One way out is to introduce 
one more pair of Bosonic fields. To perform the spatial integrations it turns out to be convenient 
to represent the Bosonic fields eft using new Bosonic fields , i.e.: 



^\x) EE £^( 



X — Xi) 



(91) 
(92) 



and using the "Lagrange multipliers" ip( a ^ (x) , x(x) , to enforce the three constraints Q92JI . At this 
point, the Gaussian variables 4>f^ are decoupled and can be integrated out. 

Skipping intermediate steps and using the fact that all the replicas are equivalent one can write 
the Green function as a correlation function: 



G( P ,z) 



1 



lim-y d d xd d y e^-") J D[^ a \ ^), X ,x) (x)^\y)A N exp(S ) 
where we have introduced the following quantities: 

So ee i J2 I d d x^ a) {x)^ a \x) + / d d x X (x)x(x) 

a=l,n 



(93) 
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-i £ Jd d xd d y^ a \x)f(x-y)^(y) 

A = J <i d xexp 3>(x) 



a=l,n ' 



x2 

a)/ 



$( X ) = +i / _ X ) X (y) _ ' V / ( g 4) 

In order to take easily the limits N,V — ► oo in (|93j) we shall resort to a grand canonical 
formulation of the disorder, introducing the partition function Z[p] = ^ N A N p N /N\ expSrj. Since 
the average number of particles N = Vp{A), in the n — > limit we have that (A) = 1 and the 
'activity' p is just the density of points, as before. Furthermore, the Gaussian integration over the 
fields xfj^ is easily performed, leading to the field theory: 

Z [p, z]=J D[^), X , x] exp (S' + Si) , (95) 

where 3 

S' Q = ~Trhxf+ [d d x X (x)x(x)-//± £ / ' d d xd d y ^ \x)f^ l \x - y)^ (y) 

a=l,n 



s i = pJ d d x (exp$(x)) 



The resolvent is related to the correlation function by ea. ()59[) . Before computing ()59|) . let us turn 
to the symmetry due to the translational invariance. Since e n (i) = constant is an eigenvector of 
zero eigenvalue of the matrix @ for every disorder realization, we see that Eq. (|15|) implies: 

G(p = 0,*) = ~. (96) 

Interestingly enough, in the framework of the field theory introduced above, that constraint is 
automatically satisfied, due to the Ward identity linked to that symmetry. The interaction term 
<£(x) is indeed invariant under the following infinitesimal transformation of order rj: 

6fo l \x) =in(z + X {x)) , (97) 
5 X (x) = 2in I d d yr\x - y)^ (y) (z + X (x)) , (98) 



while the whole variation of the action S = S' + S\ is due to the non-interacting part S' : 



5S = -"?J— I d d x^ (1) (x) (99) 



The invariance of Si hence leads to the following Ward identity: 



<* + *(*)> = j^M {X) i x ) I d d y^ 1] (y)) (ioo) 

That is enough to prove that 

Vy^O^Mfo) (101) 



diverges as 1/z at small zeta. If we combine the previous result with the exact relation (x(x)) 
—pf{0), that can be derived with a different argument, to obtain: 



z-pf(0))^ = I d d y{^(x)^(y)) 



(102) 



3 The quantity In / is computing considering / as an integral operator 
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that, together with ()59[). implies the expected constraint G(p = 0, z) = 1/z. 

The interacting term in (|96j) is complicated by the presence of an exponential interaction, mean- 
ing an infinite number of vertices. In order to perform the explicit computation of the resolvent 
G(p, z) one has to introduce some scheme of approximation. We have chosen to deal with the high 
density limit, where many particles lie inside the range of the interaction f(r). The high density 
limit {p ^ 1) of (|59|) is the typical situation one finds in many interesting physical situations, 
for example the glassy phase. In order to extract the leading term let us make the expansion 
exp$(x) — 1 ~ &(x). In that case the integration over the fields %iX ls trivial, because of: 



D[x] exp (^J d d xx(x) 



9 I d d yf(x -y) + X (x) 

l[6(x(x) + pf(Q)) (103) 



and the fields ip^ are free. In fact, introducing the quantity a = z — pf(0), one remains with the 
Gaussian integration: 

Zoc [ Df^exp-i [ d d xd d yY,^ {a) (x)K- 1 (x-y)iP < - a) (y) , (104) 

where the free propagator K is defined by: 

K-\x-y) = r\x-y) + -S{x-y) . (105) 

a 

It is then easy from (|59|) and (|104|) to obtain the result: 

G (p,z) = 3_, (106) 

z - e[p) 

where 



<P) = \Jptftp) ~ /(0) (107) 

We see that at the leading order, a plane wave with momentum p is actually an eigenstate of the 
matrix M with eigenvalue E, and the disorder does not play any relevant role. In other words, 
inside a wavelength 2tt/p there is always an infinite number of particles, ruling out the density 
fluctuations of the particles: the system reacts as an elastic medium. 

Let us finally obtain the density of states at this level of accuracy, using Ea (|18jl and Go: 

g E (E) =6(E-pf(0j) . (108) 

We obtain a single delta function at p/(0), that is somehow contradictory with our result for the 
dynamical structure factor: from the density of states one would say that the dispersion relation is 
Einstein's like, without any momentum dependence! The way out of this contradiction is of course 
that in the limit of infinite p both e(p) and pf(0) diverge. The delta function in eq. (|108|) is the 
leading term in p, while the states that contribute to the dynamical structure factor appear only in 
the subleading terms in the density of states. The same phenomenon is present in the simpler case 
discussed in the previous section. 
subsectionOne loop 

We have seen above that, with the expansion of exp<I>(x) up to first order in <&(x) the fields ip 
are non interacting and no self energy is present. Now we shall see that the one- loop correction to 
that leading term provides the - contribution to the self energy. In fact by adding the quadratic 
term to S\ the total action becomes (in the n — * limit): 
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S = S' + S 1 ~S + pJ d d x U(x) + ^$ 2 (x)) 



£ I d d xd d y X (x)f^(x-y)x(y)+i I d d xb(x) X (x) 



e(2) 



l -Y. J d d xd d y^(x)c{x-y)^ a \y) 



a=l,n 



/ d d x 



Ea=l,n(^ (a) (^)) 



2\ 2 



(109) 



where 



*(x)+p/(0)-| / d a y 



b(x) 

c(x-y) = f [ -~ 1 \x-y) + 



z + x(x) 



-f(y-x) 



6(x - y) 



(110) 



After doing some computations and adding the two contributions coming two different diagrams 
one gets 



G(p,z) 



G (p, z) + G 2 (p,z)- 



1 r d d q 



pj {2tt) c 



2 G Q (q, z) ( pf(p - q) - pf(q) 



(111) 



The Dyson resummation of all the higher orders terms, that is built by 'decorating' recursively all 
the external legs Go with the one loop correction in (|lllj) gives: 



G(p,z) 



1 



z - e(p) - Ei(p,z) 



where the self-energy Xi(p, z) is given by 

1 



£i(p,£) = 



pj (27T) 



— % G (q, z) pf(p - q) - pf(q) 



(112) 



(113) 



Let us study in details the low exchanged momentum limit of Ea. (|113|) . It is clear that at p = the 
self-energy vanishes, as required by the Ward identity ()100|) . We need to expand f(p — q) for small 
p, that due to the spherical symmetry of / yields 



fip-q) = Rq)-(p-q) — + 0(p 2 ), 

Q 

= f(q) + (p-q) € -^ + 0(p 2 ). 

qp 



(114) 
(115) 



Substituting ()115|) in (|113|) . and performing explicitly the trivial angular integrations in dimen- 
sions d we obtain 



Si (p,z) ss p* 



-,1-d 



pdir d / 2 T{d/2) ,/(, 
2 i-d 



dqq° 



P 



pd-K d l 2 T(d/2) Jo 



i W(q)f 

z - e(q) 
d-l 



q'(e)(z-e) 



(116) 
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Figure 4: The diagrams of the 1/p expansion that are taken into account in our approach. The 
numbers correspond to the particle-label repetitions 

In the last equation, we have denoted with q(e) the inverse of the function e(q). Setting now 
z = E + i0 + , and observing that e(p) ~ Ap 2 for small p, we readily obtain 

ReX 1 (p,E + iO+) « p 2 , r^de ■ (U7) 



pdTr d / 2 T{d/2) Jo q'(e)(E-e) 
I m ^E + ^) „ - - j - m - m Aq( E)f . (118) 

Since the principal part is a number of order one, the real part of the self-energy scales like p 2 
(possibly with logarithmic corrections), and thus the speed of sound of the system renormalizes due 
to the 1/p corrections. As a consequence, the function q(E) is proportional to E 1 ! 2 ~ p at the 
maximum of the function of p Se(p,E), and the width of the peak of the Se(p,E) will scale like 
p d+2 . It is then easy to check (see (|90|l ) that in frequency space the width of the spectral line will 
scale like 

r oc P d+1 , (ii9) 

as one would expect from Rayleigh scattering considerations. 

The result (|119|) for the asymptotic regime p « 1 has been found at the one loop level. In order 
to predict correctly the spectral properties at very low external momentum p, it turns out that one 
must study the behavior of the two loop contribution, that can be done in details. Nevertheless, the 
one loop result is already a good starting point to perform detailed comparisons with the numerical 
simulations. The disadvantage of this approach is that it works near band edge (uj = at high p 
but is not suited for producing the whole spectrum. Due to the complication of the action it is not 
clear how to do a variational computation and the best it can be done at the present moment is a 
CPA-like approximation described in the next section. 

6.3 A CPA like approximation 

As in the previous section we consider the resolvent G(p, z) Our aim is to compute G using the 
appropriate self-consistent equations. A partial resummation of the expansion for the resolvent can 
be written as 

z- \{p) -E(p, z) 

The self-energy S(p, z) is then expanded in powers of 1/p [21 EJ in the relevant region where p = O(z). 

If we reformulate the 1/p expansion in a diagrammatic way we can identify those diagrams with 
the simple topology of Fig. 0] Topologically, these diagrams are exactly those considered in the 
usual lattice CPA and in other self consistent approximations. The sum of this infinite subset is 
given by the solution of the integral equation: 

£(P, z) = - /7S3 [p (/(<?) " f(P ~ 9))] 2 G(p, z), (121) 

where the resolvent is given by Eq. 1)120(1 . The solution gives us the resolvent, and hence the 
dynamical structure function and density of states (Eg. 11221 below). 
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Figure 5: Top: dynamic structure factor as obtained from Eq. (jl21j) (full line) and from simula- 
tions 7 (dashes). Bottom: the density of states divided by lo 2 (Debye behavior) as obtained from 
Eq. ljl23j) (full line), simulations (dashes), and first order in the 1/p expansion (dots). 



We are interested to study the solution of Eq. ()121|) for different values of z and p. To be 
definite, we consider an explicit case where the function f(r) has a simple form, namely f(r) = 
exp[— r 2 /(2cr 2 )]. This is a reasonable first approximation for the effective interaction [H]. We shall 
take a as the unit of length and set po = 1/u, that is a reasonable choice for po for this Gaussian 
/(r), as discussed in 8 . In this particular case we will solve numerically the self-consistence equa- 
tion. We will also evaluate by simulation (using the method of moments |41j ) the exact dynamical 
structure function and the density of states by computing the resolvent for concrete realizations 
of the dynamical matrix, considering a sufficiently high number of particles so that finite volume 
effects can be neglected. These numerical results will be supplemented by analytic results, that are 
/-independent and can be obtained in the limits p — > oo and p — > 0. 

The infinite momentum limit is particularly interesting because of the remarkable result [01 E] 
that the density of states g(oS) can be written as 

9M _ lim " 2 f'(y>. (122) 

yv ' P^oo k B Tp 2 v ' 

We easily find that in this limit Eq. ()121j) can be written as: 

_L_ _ £ _ /(0) _ AQ(Z) _ jJ0-P( q)G ( q , z) , (123) 



24 



where 



Q{z) = lim G(p,z) (124) 



and 

A=(2t:)- 3 lf 2 (q)d*q. (125) 



A simple approximation consists in neglecting the last term in the r.h.s. of (|123|) . that is 
reasonable at large z. This approximation implies a the density of states that is semicircular as a 
function of a> 2 , with width proportional to ^fp and centered at uj 2 = pf(0). Translational invariance 
also requires low-frequency modes. These are given by the neglected term, and in fact it is easy to 
show that at high density it produces a Debye spectrum that extends between zero frequency and 
the semicircular part. 

In the limit p —* 0, the leading contribution to S" comes from g > p in Eq. ()121|) . where 
G(q,z) ~ G(z), so we can write for the peak width T{p) w Tq(j>), where 

g{u p ) f d 3 q 



r b) = Kp^^jr 7(^3 [/(«) - /(p - 9)J • ( 126 ) 

The integral is of order p 2 , so if the spectrum is Debye-like for small frequencies, we get T(p) ~ p 2 . 

These considerations are verified by the numerical solution of the Gaussian case, that are shown 
in Fig.EJfor p = 1.0<r -3 together with the results for the simulations 0. Note the good agreement, 
to be expected for high-densities, and how, for large p, S^\p,uj) (Fig. 03 top) tends to the density 
of states. The density of states from the self-consistent equation (Fig. ^bottom) also agrees very 
well with the results from simulations, and is a big improvement over the first term of the expansion 
in powers of The two contributions (Debye and semicircle) mentioned above can be clearly 
identified. As expected our approximation fails in reproducing the exponential decay of the density 
of states at high frequencies, that is non perturbative in 1/p [3] and corresponds to localized states. 

Next in Fig. El we plot the linewidth as a function of p as obtained from Eq. (jl21j) . Notice that 
we recover the behavior predicted from the first two non-trivial terms in the expansion in powers 
of p~ l jjj: the linewidth is proportional to p 2 at small p (also predicted by the argument above), 
then there is a faster growth and finally it approaches to a constant as S^ 1 ' (q, uj) starts to collapse 
onto the density of states. The inset shows that the contribution Eq. (|126|) is indeed dominant 
at small p. However, accurate p 2 scaling is found only for very small momenta (p/po < 0.1), 
while experiments are done at 0.1 < p/po < 1- In this crossover region, our approach predicts the 
existence of non- universal, model dependent small deviations from p 2 , that are probably hard to 
measure experimentally. In any case, the effective exponent is certainly less than 4, in contrast with 
lattice models and consistent with experimental findings. Similar conclusions can be drawn from 
mode coupling theory (see fig. 8 of 



6.4 The disappearance of phonons 

New phenomena are present the function / is no more positive or in the vectorial case. In these 
case it is possible that the spectrum arrives up to negative values and the density of states qe does 
not vanish at E = 0. In principle for a given model we should not expect a sharp transition from 
the situation where the Debye spectrum holds and gE{E) oc E 1 / 2 because tails are always present. 
However often tails are small and one can effectively observe a transition from the two regimes. 

In the framework where the density of states is computed from a simple integral equation, the 
tails are neglected and we are in the best situation to observe such a phenomenon. 

Let us call r a parameter that separate the two regimes. A detailed analysis show that a r = 
the density of states behaves as 

g E (E) oc E l '\ (127) 
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Figure 6: Peak width vs. p, for p = 0.6a 3 . The inset shows that Tq(p) is the dominant contribution 
at small p. 
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while at small r and E we have 

g e (E,r) = E 1 / 4 h(rE 1 / 2 ). (128) 

The detailed argument is a delicate but let us consider an heuristic version of it. 

With some work it is possible to obtain from (|121j) an integral equation even for the density of 
states. As a matter of fact, defining Q(z) = G(p = oo, z), the density of states turns out to be 

g E (E) = lmg(E + iO+). (129) 

Q being the solution of the following equation: 

' 1 f(0)-AG(z)-B(z) (130) 



pg{z) p 

where 



A =(2^ jf\q)d\, 



B{z)= J-0^f 2 (q)G(q,z) , (131) 

With this equation, one needs to know the resolvent at all q to obtain the density of states, due to 
the last term in the r.h.s. This can be done by solving numerically the self-consistent equation of 
ref. but here we perform an approximate analysis, that is more illuminating. 
The solution of the previous equation is 



r( , + v / «(^) 2 ~ 4Ap 

9(z) = ^ , (132) 

where 

a(z) = /(0) - - + B(z) (133) 
P 

The crudest approximation is to neglect the dependence of B{z) from z (we will also assume that B 
is a smooth function of the other parameters). In this case Ea. ll3Ul is quadratic in Q, and one easily 
finds a semicircular density of states. But the semicircular spectrum misses the Debye part, and a 
better approximation is needed. So we substitute G in the last term of the r.h.s. by the resolvent of 
the continuum elastic medium Go(z,p) = (z — £?(p)) _1 . This is reasonable because the f 2 (q) factor 
makes low momenta dominate the integral, and due to translational invariance G(z,p) ~ Go(z,p) 
in this region We shall be looking at small E, so to a good approximation 



B(z) = B + iB lZ 1/2 (134) 
We have two limiting cases, depending on the sign of of a(0) 2 — \Ap. 

• When the semicircular part of the density of states does not reach low frequencies, the square 
root can be Taylor-expanded, and one gets 

g e (E) « E 1 / 2 , (135) 

that is precisely Debye's law. 

• In the opposite situation, on the other hand, the semicircle arrives also at negative values of 
E and qe{E) is different from zero also if B% where zero. 



Exactly at the critical point we get gE(E) oc y \[E = E 1 ^ that is the announced result. 
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Mathematically, the instability arises when G(0) develops an imaginary part. This can only come 
from the square root in previous equations. Notice that this instability is a kind of phase transition, 
where the order parameter is Im<7(0). Doing a detailed computation on finds that this order 
parameter behaves as t@, with /3 = 1/2. 

Since the behavior of the propagator at high momentua does not strongly affect the dispersion 
relation we do not expect deviations from a linear dispersion relation. This has been checked 
either numerically or solving numerically the self-consistency equation given a particular choice for 
the function f(r) (see the numerical results section below). It can be argued that this kind of 
phenomenon is responsible of the Boson peak |12| . however we cannot discuss this point for lack of 
space. 

sectionCorrelated points 



6.5 Various models 

When the points are correlated things become more difficult. Already it is difficult to study the 
statistical properties of correlated points and it is more difficult to study the properties of the 
matrices that depends on these points. Although some approaches have been developed that allow 
us to deal with two points correlations [El Ej , it is not evident how to treat the general case where 
many points correlations are present. 

A particular interesting case is when the matrix and the distribution of the points are related. 
In the best of the possible words there should be extra symmetries that express this relations. 

This field is at its infancy, so that I will only describe some general results, without presenting 
applications, that for the moment do not still exist for the case of Euclidean random matrices. 

In the general case we consider an Hamiltonian i2"[a?], such that the stationary equations 

dH 

Ft[x] ^dx~ = ° (136) 

have a large number of solutions. 

Let us label these solutions with a index a. The probability P[x] of a configuration of the points 
x is assumed to be given by 

P[x] oc y~] w a 5(x — x a ) , 

a 

w a = D[M a )) eitp(-(3H[x a }) , (137) 

where the matrix M is given by 

M >* - Hk l '=" ' <138) 

This problem arises for the first time in the framework of spin glasses |44| . but its relevance to 
glasses has been stressed for the first time in [35] The function F may selects the different type of 
stationary points. 

Different interesting possibilities are: 

• D[M] = 1, i.e. all stationary points have the same weight. 

• £)[M] = sign(det[M]), i.e. all stationary points have a weight that can be 1 or -1. 

• D[M] = 1 only if all the eigenvalues are positive, otherwise it is zero (i.e. minima are selected) 
We are eventually interested to study the properties of the matrix 

ttt4 (139) 

M[x] -z, y ' 

when the points x are extracted with the previous probability. 
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6.6 A new supersymmetry 



For simplicity I will only restrict myself to the case z = where some symmetry are present when 
D[M] = sign(det[M]). 

Indeed it is evident that 

f P[x]d[x]A[x] oc f d[x]D[M}ex.p(-pH[x])detM[x]\Y[8{Fi[x})A[x] (140) 
J J . 

The last term simplify to 

J d[x]ex.p(-/3H[x])detM[x]Y[S(Fi[x])A[x], (141) 

i 

when L>[M] = sign(det(M)). For simplicity let us assume that this is the case, without discussing 
the physical motivations of this choice. 

Using usual representations we can write 



M lk = J Mx, A, V, V>] , (142) 
where d/U is a normalizes measure proportional to 

d[x]d[X]d[^}d[iJ] exp ( -0H[x] + i J2 X kF k [x] + £ M jjh ^ k I . (143) 
V k k,j J 

Here the ip are really Fermionic variables (i.e. anticommuting Grassmann variables): they have 
been introduced for representing the determinant, however they can also be used to compute the 
matrix elements of the inverse of the matrix M. 

It was rather unexpected |46[ I53j to discover that also for (3 ^ the measure d[i is invariant 
under a transformation of Fermionic character of BRST type (a supersymmetry in short), (see for 
example |52| I47j). If e is an infinitesimal Grassmann parameter, it is straightforward to verify that 
(|143|) is invariant under the following transformation, 

Sxi = eipi 5Xi = -e(3ifji 6ipi = -exi 5ipi = (144) 

This symmetry is extremely important and it is crucial to use approximation methods that do 
not break it (1H][HT)]. This has b een stressed in the spin glass case |53j in the context of the Tap 
equations 



6.7 Physical meaning of the supersymmetry 

In this section we will follow a reasoning allowing for an intuitive explanation of the physical meaning 
of the supersymmetry in terms of a particular behavior of the solutions of the stationary equations 

It may be interesting to concentrate the attention on some of the Ward identities that are 
generated by the supersymmetry. The simplest one is 

@jil>k) = i{^jXk) ■ (145) 

Apparently the equation is trivially satisfied. Indeed it is convenient to consider a slightly 
modified theory where we make the substitution 

Fj ->Fj-h . (146) 
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The final effect is to add an extra term in the exponential equal to iXjh. In other words the r.h.s 
of equation 11451 is 

d ( x k)h | hA7 , 

-^jr lh =° (147) 

The l.h.s can also be computed and one finds that eq. (|145|) becomes 

= m U=o • (148) 

a 

If we consider solutions where the det(M) / we have that 

dw n , 



dh 



|h=0 



0, (149) 



so that the previous equation seems to be always satisfied. 

However the we must be careful if most of solutions have a very small det(M). In this case if we 
first sent N to infinity and after we do the derivative with respect to h we can find that the set of 
solutions in not a continuos functions of h. Solutions may bifurcate or disappear for any arbitrary 
small variation of h and the previous relations are no more valid. 

This phenomenon has been studied in the framework of infinite range random matrices, where 
the function that plays the role of H is the free energy as function of the magnetizations (i.e. 
the TAP free energy) in spin glass type models. One finds that depending on the parameters 
there are two phase, one where the supersymmetry is exact, the other where the supersymmetry 
is spontaneously broken. This last phenomenon has been discovered last year and at the present 
moment one is trying to fully understand its consequences. 

In the framework of Euclidean random theory there are two questions that are quite relevant 
and may be the most interesting for Euclidean Random theory: 

• How to construct and to use in a practical way a formalism where the supersymmetry of the 
problem plays a crucial role? 

• How to find out if there is a phase where supersymmetry is spontaneously broken and which 
are the physical effects of such a breaking. 

It is quite likely the response to these questions would be very important. 
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